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TRANSITION TO TURBULENCE IN PLANE CHANNEL FLOW 

By 

Sedat Biringen 
SUMMARY 

This work, performed under Grant No. NAG-1 -228 from NASA/Langley Re- 
search Center, involves a nunerical simulation of the final stages of 
transition to turbulence in plane channel flow. Three-dimensional, incom- 
pressible Navier-Stokes equations are nunerical ly integrated to obtain the 
time- evolution of two- and three-dimensional finite-amplitude disturbances . 
Computations are performed on the CYBER-203 vector processor for a 32x51*32 
grid. Results are presented for no-slip boundary conditions at the solid 
walls as well as for periodic suction-blowing to simul ate active control of 
transition by mass transfer. Solutions indicate that the method is capable 
of simulating the complex character of vorticity dynamics during the various 
stages of transition and final breakdown. In particular, evidence points to 
the formation of a A-shape vortex and the subsequent system of horseshoe 
vortices inclined to the main flow direction as the main elements of trans- 
ition. Calculations involving periodic suction-blowing indicate that 
interference with a wave of suitable phase and amplitude reduces the distur- 
bance growth rates. 



1. INTRODUCTION 

Recent experiments by Nishioka, Asia & li da (1981) have shown that 
transition to turbulence in a plane channel flow follows a sequence of 
events similar to that observed by Klebanoff, Tidstrom & Sargent (1962) in 
the bound ary- layer transition. In this work, a direct numerical integration 
of the Navier-Stokes equations is performed in an attenpt to simulate these 
events in plane channel flow* during the later stages of transition. 

In their experiments, Nishioka et al . (1981) measured the streamwise 
mean and fluctuating velocities, U x and Ul ', respectively, at a fixed 
streamwise location at a subcritical (linearly stable) Reynolds nuriber, Re * 
5000 and simulated the various stages of transition by varying the distur- 
bance amplitude. Their observations show that subcritical instability takes 
place at a threshold amplitude of (uj^/Uo = 0.01, where U 0 is the mean 
velocity at the channel centerline, the evolution of this instability is 
evidenced by the intensification of the spanwise variation of the wavefront 
which develops into a peak-valley structure. Nishioka et al. (1981) 
observed that flow development follows a trend which is similar to trans- 
ition in the boundary-layer (Klebanoff et al. 1962, Kovasznay, Komoda & 
Vasudeva 1962): local shear layers are formed away from the wall at span- 

wise peak positions (u 1 /U 0 = 0.06), these shear layers start to exhibit 
a "kink" which is the manifestation of secondary instability and is accom- 
panied by a "spike" in the uj/U,, = 0.11). in rapid succession, two-, 
three-, five-, and multi-spike stages are observed with increasing ampli- 
tude of the primary disturbance. Nishioka et al. (1981) presented evidence 
that in the final stages Of transition, the flow Starts 1 to develop struc- 
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tures very similar to those found in fully developed wall turbulence. 

During this stage, the flow field is characteri zed by the development of a 
viscous sublayer, occurrence of the typical "streaks" close to the wall and 
the formation of horseshoe vortices sometimes referred to as the building 
blocks of wall turbulence (Theodorsen 1954). The present work simulates 
this sequence of events. 

Effects of three dimensionality on transition have first been document- 
ed in detail by Klebanoff et al. (1962). Accordingly, three-dimensionality 
manifests itself mainly in the spanwise velocity variations resulting in the 
production of streamwise vorticity which, in turn, interacts with the span- 
wise vorticity and drives the flow to breakdown. Orszag & Kells (1980) and 
Patera & Orszag (1981) have expanded on this idea and studied the suscepti- 
bility of plane channel flow to three-dimensional disturbances by numerical- 
ly integrating the three-dimensional Navier-Stokes equations. Their compu- 
tations at subcritical Reynolds numbers revealed some interesting aspects of 
subcritical transition. They found that initially two-dimensional distur- 
bances which are finite-amplitude two-dimensional Orr-Sommerfeld eigen- 
solutions, decay slowly whereas the interaction of two-dimensional and 
three-dimensional disturbances drives the flow to instability for Reynolds 
numbers as low as 1000. In the present work we focus on the three-dimen- 
sional and nonlinear mechanisms that characterize transition in the final 
stages and during flow breakdown. Hence, we employ the full three-dimen- 
sional time-dependent Navier-Stokes equations. Calculations were performed 
at a linearly stable Reynolds number (Re = 1500), with finite-amplitude two- 
and three-dimensional eigensolutions of the Orr-Somerfeld equation used as 
the initial conditions. No attempt was made. herein to study the effect of 
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different initial conditions or of Reynolds nunbers; this is the subject of 
another investigation. We also investigated the idea of using periodic suc- 
tion-blowing at the solid walls as a means of controlling transition. These 
computations were performed at Re * 7500 which is the linearly unstable 
range for this flow. Since transition control must be applied prior to the 
flow breakdown (e.g. before the occurrence of multi-spikes), disturbance 
amplitudes during these calculations were not large enough to impose stiff 
restrictions on the allowable time step. 

In section 2, of this report, the nimerical methods used in this study 

are briefly discussed. In section 3, results of calculations are presented 

and finally in section 4 a summary of results and some concluding remarks 
are given. 




2. THE CALCULATION PROCEDURE 

The calculation procedure is based on the incompressible Navier-Stokes 
equations in primitive-variable form, 


3u i 9u i 1 an* 3 u i 

_L + u. _I = - 1 1£_ + v 2_ 

at ax„ p 3x . ax # ax fl 

i i i i 


and the continuity equation, 
3u, 


= 0 


ax 


( 1 ) 

(2) 


where u.. are the velocities along the x.. directions, p is the density 
v is the kinematic viscosity and p' is the total hydrostatic pressure. 

The equations are non-dimensionali zed by the mean centerline velocity Uo 
and the channel half-width, h. The flow is assumed to be driven by a con- 
stant mean pressure gradient 2/Re, where Re is the Reynolds number 
given by U 0 h/v. Also, the convective terms are written in a form which 
prevents occurrence of nonlinear instability in the numerical solution 
procedure by ensuring conservation of momentum and energy (Mansour, Fer- 
ziger, & Reynolds 1978). The final form of the Navier-Stokes equations 
reads. 


au. 

To.. 

(!!i -!!i 

+ U £ 


8x i 

l 8x t 3x ( 


+J_ _!!?!_ 

ax. Re Re ax^ax^ 


(3) 


where P = p'/p +— u^u^ is the pressure head and <5^ is the Kronecker 


delta. 



The flow is assumed to be periodic in the streamwise x L and the span- 
wise X 3 directions along which the flow field variables can be expanded in 
terms of Fourier series. This enables the use of the pseudo- spectral method 
(Orszag 1972) to calculate the spatial derivatives along xi and X 3 by 
use of discrete Fourier transforms. Considering transforms in the direc- 
tion, along which there are N x equal ly spaced mesh points, the velocity 
component u 1 can be written as 

N/2-1 

ui(xi) = z ui(ki)e 1X1 (4) 

N 

n i ■ - - 
2 

where x^ - max^ , m - 0,1,..., N-l and k^ = 2 irni /N^axi . Accordingly, the 
Fourier transform of u x is 

1 Ni-1 .. 

u i ( k i ) = - E Ui (xi )e~ 1klXl ( 5 ) 

Ni m=0 

The spatial derivative of ui along xi can now be written as 


3 Ui (Xl ) 
axi 


Ni/2-1 ^ ik 

E i ki ui (ki )e 11 

ni= -Ni/2 



The derivative can be computed by forming the Fourier transform of 14 (xi), 
multiplying the result by ik 1 and computing the inverse transform. For 
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periodic functions, the pseudo- spectral method provides a means by which the 
spatial derivatives are evaluated with maximum accuracy for a given number 
of grid points. Along the x 2 -direction, a mesh stretching that concen- 
trates grid points close to the solid walls is employed (Moin, Reynolds & 
Firziger 1978). The resulting mesh enables the resolution of the sublayer 
that is formed during transition for y+ < 2 , where y+ is the coordinate 
along x 2 , inwall units. Spatial derivatives along x 2 are evaluated by 
a second-order finite-difference scheme on this non-uniform mesh. 

The governing equations were numerically integrated by the semi-implic- 
it method of Moin et al. (1978). This procedure employs the explicit Adams- 
Bashforth method for the convective terms and the implicit Crank-Nicholson 
method for pressure and for the viscous diffusion terms. In order to start 
the two time-level Ad ams-Bashforth method, the Eul er- implicit method is used 
at the first time step. 

Once the governing equations are discretized in time, a two-dimensional 
Fourier transform along the periodic directions xi and X 3 transforms the 
equations into the k 1 -k 3 wave-nunber space. The transformed equations are 
written below in block-tridiagonal form for inversion along X 2 


A f"! 1 B ff 

J+l — ~J 


+ C F ^, 1 = R y 

J ' J 



In (7), A, B, and £ are coefficient matrices, F is the solution vector 

j 

at the advanced time level, n+ 1 , and at the X 2 -directional node, j; Bj is 
right-hand side vector that contains the convective, diffusive and pressure 
terms at the previous time levels. These are given as 
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m m 

w 

*Ji 

U2 

xv 

A * 

C2 J 

0 

0 

C2 J 

0 

0 

1 - 

0 

0 

U3 

XV 

p 


0 0 

--- » 

0 

0 

C1 J 

C2 d 

0 

-Re. Cl. 

J J 


B = 


C = 


B 2 j-* 

0 

0 

±(ki)jRe 

0 

B2.-* 

0 

± (k 3 ) -Re 


*( k 3)j 

81 j 

0 

0 

0 

0 

B2 .-<f> 

-Re.Bl . 



J 

J 

A 2 . 
J 

0 0 

0 


0 

A 2 . 0 

0 J A 1 

0 


0 

0 


0 

0 A 2 J 

J 

-Re.Al . 
J 



and 


R, = 0 


R 2 = Re | - 2u * + / 92u i + 3 2u i Vn 
aT V*i 2 3 x 3 2 


n 


+ RelL 

3X1 


- 2R el I Hj - I Hi"”^ 


3 2 ui' 


3X 2 ‘ 


R 3 = Re I - 2u * + / 92 U 2 + 3 2u 2 ^ n 


at 


3 Xi 


9 X 3 ‘ 


n 


♦ Re iL 

3 X 2 


- 2 Re [1 H'J - A hJ 

l 2 2 


1 u n ’M 92 Mi 


3 x 2 ' 


R4 = Re 


2u? +/3 2 u 3 + 32 u 3 Xn 

AT \ 9 xi 2 3 x 3 2 


n 


+ Re^iL 

*3X3 


. 2Re[l H 3 n - 1 H 3 1 " 1 


32 u? 

3 X 2 2 
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(no summation over i) 


also, 


♦ s 


2Re +(4 ♦ 4' 

AT 1 


H 1 ■ - V 


' 8U J_ _^v 

<3X. 3x i , 


and, 


R. = R! 
J J 


Coefficients of the finite-difference operators that appear in the matrices 
A, and C_ are given as 


C2 j ■ 2/ Vi<Vi ♦ "j)- B2 j ■ - 2/ ViV ■ 2/ Wi ♦ V 

C1 j * - A1 j ■ - v 

A j+1 = ( X2 )j+1 " ( X2 )j 

Since all the flow variables in the solution vector contain an imaginary and 
a real part, the block-inversion process is applied twice for each pair of 

k x and k3, which are the wavenumbers along xi and X3, respectively. 

The assumption of periodicity in x x and x 3 eliminates the necessity 
of applying explicit boundary conditions along these directions. However, 
due to the presence of solid boundaries along the x 2 direction, no-slip 
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boundary conditions are imposed on u x , u 2 , and u 3 and the pressure at 
the wall is calculated by a second order approximation from the interior of 
the flow field. That the pressure boundary conditions are consistent with 
the x 2 -momentum equation at the wall, has been shown by Moin et al. (1978). 

Initial conditions were prescribed from the two- and three-dimensional 
eigensolutions of the Orr-Sommerfeld equation by considering that even for 
subcritical Reynolds numbers, plane channel flow can be driven to instabili- 
ty if the least stable two-dimensional finite-amplitude Orr-Sommerfeld 
eigenmodes are interacted with finite-amp! itude three-dimensional eigenmodes 
(Orszag & Kells 1981). The most explosive situation arises when the three- 
dimensional eigenmodes are aligned with the main flow direction at ±45 to 
±60 degrees. Accordingly, we have used the following initial condition: 

±j(x) = U(x 2 ,0,0) + u 2D (x 2 )e iax i + u 3D (x 2 )e ictx i ±i3x B (8) 

Here, U(x 2 ,0,0) is the parabolic velocity profile of plane channel flow. 
The eigenfunctions u 2D (x 2 ) and u 3Q (x 2 ) correspond to two-dimensional and 
three-dimensional solutions of the Orr-Sommerfeld equation at Re * 1500, 
respectively. The two-dimensional solution was obtained for « = 1 whereas, 
the three-dimensional solution was obtained for a = 1, e = ±1. A computer 
program, which essentially uses the Kaplan filtering technique was used for 
the solution of the Orr-Sommerfeld equation (Reynolds 1967). The final 
amplitudes were chosen so that the maximum value of the x r directiona1 two- 
dimensional disturbance was set equal to 0.11U 0 and the maximum amplitudes 

of the xrdirectional three-dimensional disturbances were each set equal to 
O.O5U 0 . 
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3. RESULTS AND DISCUSSION 


The finite-difference system (7) was solved on the CYBER-203 vector 
processor at NASA/Langley Research Center. A 32x51x32 mesh was employed 
along the X!-, X 2 -, and X 3 -directions, respectively. The computer code was 
fully vectorized and vectorized library subroutines were used for the main 
computational operations that the solution technique employs. These vector 
operations are mainly one-dimensional fast Fourier transform (FFT) to cal- 
culate spatial derivatives with the pseudo-spectral method, two-dimensional 
FFT to transform the equations into k 1 -k 3 wave-nunber space and block- 
tridi agonal matrix inversion along X 2 . For the FFT operations, typical 
vector lengths were around 1000 , which is an optimal vector length to take 
full advantage of the vector processor. For the block-tridagonal matrix 
inversion (which essentially is a scalar operation), a vectorized subroutine 
that inverts a large number of tridiagonal systems simultaneously was used. 

This procedure decreases CPU time significantly by reducing the number of 
scalar operations required to invert each system separately. The fully 
vectorized code takes about 10 sec. of CPU time per time step for the 32x51x32 
mesh to solve the finite-difference system (7) on a computational box, in which 
the flow is confined between rigid walls at x 2 = ±1. Periodicity lengths 
(box lengths) along xi and X 3 were chosen so that the smallest wave 
numbers allowed in the computational domain were equal to a = 1 and 3 = 1 , 
respectively, i.e., the box length was set equal to 2 ir along these directions. 

It should be recalled that the time- advancement scheme employed in this 
work is partly explicit (on the convective terms) and partly implicit (on 
the diffusion and pressure terms). Although. in view of linear stability 
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analysis, implicit methods are unconditionally stable (extrapolation to 
nonlinear equations is sometime vague), the mixed nature df the present 
scheme as well as the time-accurate nature of the problem under investiga- 
tion necessitate adherence to stability bounds of explicit Schenes. There- 
fore, in all the calculations reported herein, the convective stability 
condition (the Courant-Fri edrichs-Lewy condition) that requires the Courant 
nunber (C.N.) to be always less than one and the diffusive stability condi- 
tion were obeyed. With (Ax 2 ) min = 0.0092 and aj = 0.025, where a T is the 
nondimensional time-step, through the course of the calculations C.N. varied 


as 


C.N. = AT 


ui_ 

AXl 


d2 

AX2 


U3_ 

AX3 


< 0.2 


(9) 


/max 


whereas the diffusive stability criterion, D, varied as 


1 

Re 


AT 


(Ax) 2 


min 


< 0.04 


( 10 ) 


so that the diffusive stability condition which requires 0 < 6.5, was also 
always satisfied. The computer program was tested by calculating the growth 
rates of smal 1-ampl itude Orr-Sommerfeld waves. For a wide range of Reynolds 
nunbers (between 1000 and 10,000) the agreement between the computed re- 
sults and the linear theory was better than 0.5%. 

In the subsequent parts of this section, results obtained from the 


nimerical integration of the finite-difference system (7) for the time-evo- 
lution of the initial disturbances are compared with the experiments of 
Ni shi oka et al . (1981) . It should be noted that this experiment was done at 
Re = 5000 (subcritical) whereas the computation was done at a lower 
aubcritical Reynolds nunber, Re = 1500. The selection of a higher Reynolds 



number (e.g., in the linearly unstable range) makes the governing equations 
very stiff and requires the use of extremely small time steps for numerical 
stability. Hence, the selection of Re a 1500 was mainly to force the 
computations into transition and breakdown with the least amount of computer 
expense. It should be noted that wall phenomena characteristic of the final 
stages of transition are essentially independent of Reynolds number as 
observed by Nishioka et al. (1981) for channel flow and by Smith & Metzler 
(1982) for boundary-layer transition. Therefore, the difference between the 
Reyn.olds numbers of the experiment and the simulation should not have any 
important consequences for the qualitative comparisons between the two sets 
of results. 

3a. Plane Channel Flow with No Mass Transfer 

In their experiments, Nishioka et al . (1981) identified the various 
stages of transition according to the number of spikes appearing in the os- 

I 

cilloscope traces of the disturbance velocity, u^ Since the time-axis of 
the experiment is interchangeable with the x^-axis of the computation, we 

I 

obtain similar traces by plotting u x along x x (over two periods) as shown 
in figure 1. In this figure, the first frame shows the sinusoidal variation 
of the initial conditions at T = 0. This is followed by the nonlineaer dis- 

I 

tortions of the initial conditions resulting in variations of u x which 
strongly resemble the oscilloscope traces at the one-, three-, and five- 

I ■" 

spike stages of the laboratory flow. In particular, the variations of ui 
with x x at T = 44 are very similar to the ensemble- averaged waveforms 
presented by Nishioka et al. (1981) at the five-spike stage. The last frame 

I 

in figure 1 shows variations of ^ at a time (T = 79) much later than the 
occurrence of the five-spike stage (T = 44). . Hence, we note the absence of 
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flotations of rich frequency content that are characteristic of fully- 
developed turbulence; this lends support to the idea that the last stage of 
transition may not be spontaneous and explosive a phenomenon as generally 

supposed. However, at this stage, the ui variations are very similar to 
velocity oscillations observed in wave packets and indicate the occurrence 
of "patches" of turbulent fluid. For the remainder of this discussion, 
results from the simulation at the spike stages (shown in figure 1) will be 
used, whenever possible, for comparisons with the corresponding spike stages 
of the experiments of Nishioka et al. (1981). 

In figure 2, we show a history of the time- evolution of the flow in 
terms of the maximum ampl itude of the two dimensional primary disturbance, 
its two-dimensional harmonic and the three-dimensional primary disturbance. 
The trends displayed by these quantities are generally similar to the re- 
sults of Orszag & Kells (1980) which they obtained from computations per- 
formed at Re = 1250. The main features of these trends are the rapid de- 
crease in the two-dimensional primary-wave amplitude and the rapid increase 
of the amplitude of its harmonic. Also, the three-dimensional primary-wave 
amplitude first increases, then gradually decreases at around T = 30. In 
the present calculations, we observe that variations of the amplitudes start 
to fluctuate as early as T = 15 but, even at later times, the fluctua- 
tions do not display an explosive trend. That no "explosive" instabilities 
were found in the present computations is in accord with the findings of 
Nishioka, As ai & Iida (1980), which imply that breakdown in channel flow is 
as gradual as the growth of instabilities found in free shear flows. 

In figure 3, we present plots of maximum root-mean- square (rms) anpli- 



tudes of ui , (ui ms ) max , over two periods along X3 . Here, we define the 
rms value as an average over Xj^ . We note that at 1*0, the wave pat- 
tern is sinusoidal and corresponds to peaks at X3 = 0, 2*, 4* and to 
valleys at x 3 = ir» 3ir. Subsequently, at later times the nonlinear dis- 
tortions of the wave front result in minima occurring at the peaks and 
maxima occurring at the valleys in accordance with the experiments of 
Nishioka et al . (1980). At later times, an increase in the frequency of the 
peak-valley structures is clearly evident suggesting an increase in the nun- 
ber of characteristic vortex structures along X3. 

Plots of velocity profiles averaged over the x 1 -x 3 plane, <Ui>, are 
shown in figure 4 for laminar (initial), late transition and "early turbu- 
lence" stages at T = 0, T = 44 and T = 79, respectively. The <ui> 
distribution at T = 79 has a strong resemblance to the turbulent 
channel flow profile, with increased velocity gradient at the wall and with 

a full profile indicative of turbulent mixing. Although, as expected, <ui> 
profiles do not show any fluctuations (or inflections), plots of instantan- 
eous velocity profiles do show very strong inflections, especially in the 
regions close to the walls (figure 5). This indicates that the interac- 
tion of two- and three-dimensional waves close to the walls is the central 
mechanism that drives the flow to instability. This is in accord with the 
idea that the flow will undergo transition only for a selected band of span- 
wise wave numbers, the most "dangerous" of which result in three-dimensional 
disturbances with maxima occurring close to the walls (Orszag & Kells 1980). 
In figure 6, the velocity profile <u+> = <ui>/u T versus y+ * X2 u T /\> is 

plotted; here u T is the friction velocity and is calculated from d<ui>/dx2 
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at the wal 1 times 1/Re. Although, the plots indicate the formation of a 
sublayer, and the change from T = 44 to T = 79 shows a gradual approach 
to the 1 aw-of-the-wal 1, the difference is still especially apparent in the 

i 

logarithmic region. At T = 44, the Reynolds nunber based on friction 
velocity is equal to 69 and, as expected, is larger than its initial 
(1 aminar) value of 55. 

Plots of pi ane- averaged fluctuations intensities, <(u 1 - <u 1 >)2>, 
are shown in figure 7 at various T. There are several interesting features 
of this figure. Firstly, at T = 17, the shift in the position of peak 
amplitude towards the channel center (x 2 > -0.6), as well as the increase 
in the maximum amplitude indicate that the development of the computed flow 
field is compatible with experimental observations pertaining to the one- 
spike stage of the transition process (Tani 1969). It will later be shown 
that during this stage there is a substantial increase in the span Wise 
vorticity, u> 2 , away from the lower wal 1 around x 2 = - 0.6. Secondly, in 
accordance with the laboratory flow of Nishioka et al. (1981) at later 
stages, the computed intensity profile displays a second peak occurring 
close to the wall associated with turbulence production. At T = 44 cor- 
responding to the five-spike stag!, we see that peak intensity has reached a 
value typical of turbulent channel flow; however, the peak occurs uncharac- 
teristically away from the wall. Finally at T = 79, the peak in the .inten- 
sity profile has moved towards the wall but even at this stage, the distri- 
bution has not ass uned the asymptotic f ul ly-turbulent form. In figure 8, 
plots of plane-averaged shear stress, Cuj - <u 1 >)u 2 >, profiles are shown. 
The increase in the magnitude of shear stress from T = 0 to T = 17 clearly 
indicates the effects of nonlinearity in transferring energy from the mean 


16 



flow to the fluctuating motion. At the five-spike stage (T = 44), the 
maximum shear stress has attained a value typical of turbulent channel flow. 
However, the locatiion of the maximum is away from the wall and roughly cor- 
responds to the location of the peak in the respective intensity profile 
(figure 7). It can thus be inferred that at this stage of the comutation, 
the energy-exchange mechanisms of fully-developed wall turbulence is not yet 
reflected by the pi ane- averaged velocity correlations. 

Spanwise variations of are plotted at various distances along x 2 
in figures 9a-9e. In figure 9a, we show initial distributions which are 
highly three-dimensional. At T = 17 (figure 9b), the distributions indicate 
a velocity defect both at peak and valley as well as a velocity excess in 
between. This is indicative of intensification of the initial streamwise 
vortices and appearance of weaker streamwise vortex pairs in accordance with 
the Benney-Lin (1960) theory. Subsequently, figures 9c-9e display the form- 
ation of additional vortex pairs indicating transport of energy down the 
wave-nunber spectrum. The spanwise synmetry imposed by the initial condi- 
tions is retained through the f ive-spi ke stage. At this stage, an estimate 
of the flow field resolution can be obtained from the spanwise distance be- 
tween peak positions, x z . Nondimensional ised by u 2 and ^ typically x z = 
115. This is larger than but comparable to X = 80, which is the typical 
spanwise length in the laboratory flow during the five-spike stage (Nishioka 
et al. 1981). It should be noted that the spanwise characteristic length in 
wall turbulence is about 100. Therefore, it could be asserted that at this 
stage present results are representative of initial wall turbulence. 

A more detailed description of the transition process can be obtained 
from contour plots of equi-shear lines, 3 u ly (ax 2 (which correspond to 


17 



approximate spanwise vorticity, <^) in the x x - x 2 plane at the position of 
maximum u^Ug. Results from the computation that correspond to the various 
stages of the laboratory flow are presented in figures 10-15 between the 
lower wall (x 2 = -1.0) and channel center (x 2 = 0.0). In figures Ua 
13a, figures 4-6 of Nishioka et al. (1981) are also shown for a qualitative 
comparison with the present results. In figure 10, contour plots corre- 
sponding to the initial conditions and in figure 11, contour plots corre- 
sponding to the "one spike" stage (T = 0 and T = 17, respectively) are , 
shown. In both the laboratory flow arid the computation, the typical head, of 
the shear layer appears very clearly at the one-spike stage, indicating the 
formation of a shear layer away from the wall at about x 2 = - o.tj due to the 
induced velocity from the streamwise vortex system. In addition, the sud- 
den dip of the shear layer from the high-velocity outer flow to the lovf- 
yelocity region clearly appears as a kink in both the computation and the 
laboratory flow. Since the grid points are finely clustered along x 2 
close to the wall, vorticity concentrations in this region are also ade- 
quately resolved by the numerical simulation. 

Figure 12 shows equi-shear lines at T = 27 corresponding to the three- 
spike stage of the laboratory flow. Due t° the secondary instability mani- 
fested in the previous stage, breakdown of flow structures into smaller 
scales is observed in both the experiment and the computation. The growth 
of the kinked portions of the equi-shear lines into the so-called "hairpin 
eddies" is clearly depicted in the experiment. The computation displays a 
similar evolution: the head of the shear layer is lifted up towards the 
channel centerline, and the kink in the shear layer i$ quite apparent in 
this stage of the simulation. Simultaneously with this activity taking 
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place in the outer (high speed) portions of the flow field, both the 
experiment and the computation show an intense shear layer developing close 
to the wall, which is indicative of turbulence generation. It is generally 
agreed that hairpin eddies which are lifted towards the centerline erupt 
into turbulent spots. However, figure 12 indicates that wall turbulence 
may also be closely associated with vorticity dynamics simultaneously taking 
place with the eruption of hairpin eddies. Contours of equi-shear lines at 
T = 44 corresponding to the five-spike stage of the experiment are shown in 
figure 13. In both the experiment and the computation, the intense shear 
layer developed in the wall region is closely discernible. It should also 
be noted that, the spanwise position of maximum u> z which appears very 
close to the wall and does not necessarily coincide with the spanwise 
location of maximum Uj/Ug where the contours are presented. The most sig- 
nificant feature of figure 13 is the existence of distinct vortex structures 
in the wall region both in the laboratory flow and in the computation. 

These vortices are inclined to the main flow direction at an angle that 
varies between approximately 14° to 40° and show a close resemblence to the 
energetic horseshoe vortices characteristic of initial wall turbulence (Hama 
& Nutant 1963, Klebanoff et al. 1962). It is these vortices roughly aligned 
along the direction of maximun extensional stress that are mainly responsi- 
ble for extracting energy from the mean shear (Tennekes & Lumley 1972). 

In figure 14, equi-shear lines at spanwise locations x 3 = (x 3 ) 0 + 2ax 3 
and x 3 = (x 3 ) 0 + 4ax 3 at T = 44 are compared. Here, we define (x 3 ) 0 as 
the spanwise position corresponding to figure 13 and ax 3 is the mesh size 
along x 3 . Figure 14 suggests that the shear layer is formed from a system 
of horseshoe vortices in succession such that as the first-born vortex 
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erupts into the channel center, a new one forms at the wall. As evidenced 
in Figure 5, the vortex lift-up is still clearly discernible at the "early 
turbulence" stage (T = 79). Here, we note the vortex lift-up towards the 
channel center but its extension is considerably shorter than that observed 
previously at the five-spike stage. This variance suggests the existence of 
different vortex structures characteristic of these two stages. In fact, we 
shall later show that the typical horseshoe vortex of 1 ate transition is not 
the basic structure of early turbulence. 

As we have noted before, the mechanism that is responsible for the 
generation of vorticity concentrations is usually explained as vortex 
stretching (and deformation) by the mean flow. The stretching and deformed 
layer moves downstream with a translation velocity that induces lower local 
velocity of the upstream edge of the vorticity layer than its downstream 
edge (Komoda 1967). An examination of Ui contours of (figures 16 and 17) 
along with the approximate spanwise vorticity contours at T = 27 (figure 

12b) and at T = 44 (figure 13b) in the Xl - x 2 plane clearly shows a similar 
trend. We observe that the nose of the vorticity layer is generally associ- 
ated with higher velocities than the upstream region and large variations of 
the local velocity exist within the layer. 

Finally, normal velocity, i^, contours are shown in figures 18-21 
during the three-spike, multi-spike and early-turbulence stages at the 
position of maximum ^/Uq. Figure 26 shows u 2 -contours in the x x - X 2 
plane at T = 27, corresponding to the three-spike stage of the laboratory 
flow. In this figure, there is clear evidence of the beginning of an alter- 
nating up-and-down flow similar to the pattern described by Kovaszny et al. 
(1967), which is characterized by the intense updrift accompanied by fluid 
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drifting down at both sides. This alternating structure of up-and-down flow 

In the xi - x 2 plane is more evident at T =44 (figure 19), corresponding to 
the multi -spike stage of the experiment. Also in this figure, alternating 
regions of fluid with scales smaller than those at T = 27 are depicted. In 
figure 20 and 21, normal velocity contours are shown in the x 2 - x 3 plane at 
T = 44 and T = 79 corresponding to the multi-spike and early-turbulence 
stages, respectively. These contours are plotted in the region between the 
lower wall (x 2 = -1.0) and x 2 = -0.92, display the evolution of alternating 
structures very similar to the characteristic streak-like structures found 
in the wall region of turbulent channel flow (Moin & Kin 1982). 

3b. Active Control of Transition by Periodic Suction-Blowing 

In this section we demonstrate the usefulness of employing periodic 
suction-blowing boundary conditions at the solid walls for controlling 
transition. All these calculations were performed at a Reynolds number well 
within the linearly unstable range, i.e. at Re = 7500. The main idea here 
is to cancel or modify the "most dangerous" disturbance by interfering with 
a wave of suitable amplitude and phase. This idea has been proven to be 
/effective in controlling boundary-1 ayer transition in an experiment in which 
the two-dimensional Tollmien-Schlichting wave was effectively cancelled by a 
control wave generated downstream of the initial disturbance. The amplitude 
and phase of the control wave (produced by a second wave generator pi aced 
downstream of the first) was adjusted to minimize the maximum disturbance 
from the first wave generator (Milling 1981). Expanding on this idea and 
assuming the two-dimensional primary wave to be the most dangerous disturb- 
ance, in this work we prescribe the boundary conditions at the solid walls 
such that 
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(a) the resulting control wave is compatible with flow periodicity 
along and Xg so that there will be no net mass transfer into 
the flow field; 

(b) the resulting control wave will be 180' out of phase with the two- 
dimensional primary wave; 

(c) the control wave amplitude is adjusted for maximum attenuation of 
the primary disturbance. 

In figure 22, results for a case in which the initial disturbances are 
two-dimensional solutions of the Qrr-Sommerfeld equation are shown. The 
maximum amplitude of u^ was set equal to 1 % of the channel centerline 
velocity. Periodic suction-blowing boundary conditions were prescribed at 
each wall at T = 5.3 for one time stop, before and after which the channel 
walls were made impermeable via the no-slip boundary conditions. The most 
effective wave produced by the suction-blowing boundary conditions was found 
to be 180 out of phase with the two-dimensional primary u 2 (corresponding 
to u 2 (1,0) in the Fourier space) with twice its amplitude. The effect of 
the control wave to reduce the growth rates of both (u|) ^-primary and its 
harmonic is evident from figure 22. In fact, the control wave has resulted 
in a decaying harmonic preventing energy transfer down the wave-nunber 
spectrum, consequently inhibiting the proper development of the transition 
process. 

To test the effects of three-dimensional ity a similar calculation was 
performed in which three-dimensional solutions of the Orr-Sommerfeld equa- 
tion were also included into the initial disturbance field. For this ease, 
the max imun amplitudes of the two- and three-dimensional disturbances were 



set at 3% and 1% of the channel centerline velocity, respectively. The most 

I 

effective control wave was found to be 180° out of phase with ( u 2 ) 20 
2.5 times its amplitude. The characteristics of the resulting flow field 
are shown in figure 23, after the periodic suction-blowing boundary 
conditions are applied at T = 5.3 for one time step. In summary, the 

I I 

control wave attenuated the growth rates of ( u i) tota i and ( u i) 2D” P r " iniar 'y> 

I 

and resulted in a decaying (u^q - harmonic as in the previous case. Note 

I 

that ( u i ) 3 q ” Primary remained unchanged. 

It is apparent from the preceding discussion that the application of 
periodic suction-blowing boundary conditions on the solid walls of the 
channel has a considerable effect in reducing the growth rates of two- and 
three-dimensional disturbance fields. In practice the technique can be im- 
plemented via a suction-blowing strip along the body span. The spectrum of 
the disturbance field just upstream of the strip must be available in order 
to activate the correct boundary conditions and obtain the most effective 
control wave. Further work is necessary to evaluate the frequency-response 
and feedback delay of such a mechanism. 
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4. SIMMARY AND CONCLUDING REMARKS 
In this study, final stages of transition to turbulence in plane 
channel flow have been simulated by a direct numerical solution of the 
Navier-Stokes equations. It is found that, in spite of the limited resolu- 
tion of the 32x51x32 grid employed in the computations, the simulation is 
capable of reproducing most of the essential features of wall phenomena 
observed in the laboratory. Grid resolution in the xi and xs directions, 

along which the flow is periodic, is found to be adequate to capture the 

sequence of events that lead to early turbulence. Vorticity contours in the 

vicinity of the lower wall indicate formation of a system of horseshoe 

vortices with legs or extensions in the x x - x 3 plane comprised of 
counter-rotating streamwise vortex pairs. Our findings are also in accord 
with the Benney-Lin (1960) theory, e.g., the present computations clearly 
depict frequency -doubling of streamwise vorticity. It should be noted that 
the initial conditions used in this work are of the Benney-Lin type and 
consist of a vorticity field with a strong spanwise component and weak 
streamwise and transverse components. Therefore, the question of the origin 
of stream-wise vorticity and its precise relation to distortions of spanwise 
vorticity remain to be addressed in future work. 

At later stages of the computation, transverse velocity contours indi- 
cate the formation of streak-like structures alternating in the spanwise 
direction. Typically, the spanwise characteristic length of these vortices, 
inferred from the spanwise variations of u a , was found to be approxi amtely 

x z this is close t0 x 2 = 100 of fully -developed wall turbulence. 

It was found that during the later stages of transition, flow field 
statistics indicate the formation of a laminar sublayer; however, the 



development of the logarithmic region and consequently the approach to 
fully-developed turbulence is slow. This gradual approach to steady state 
is also reflected in the profiles of pi ane- averaged intensity and shear 
stress. 

The main deficiency of this study stems inevitably from limited spatial 
resolution and manifests itself in several ways. Firstly, at later stages 
of the computation, insufficient mesh resolution results in lower gradients 
of the mean velocity in the viscous sublayer. This, in turn, causes less 
turbulence production in the wall region as evidenced by the absence of 
peaks in the intensity and shear-stress profiles close to the wall. Second- 
ly, the finite cut-off wave number's along x x and x 3 preevent the formation 
of a proper wave-number spectrum. As a consequence, pollution of Fourier 
modes and accumul ation of excess energy at low wave numbers become very 
significant sources of accuracy at large T. Hence, proper simulation of 
transition beyond the early-turbul ence stage necessitates the use of higher 
grid resolution. Even then, the incorporation of a mechanism to account for 
the subgrid scale turbulence would be required for a more accurate and 
realistic representation of the flow field phenomena. 

This work was supported by NASA/Langley Research Center under Grant 
No. NAG-1-228. The author is indebted to P. J. Bobbitt and W. D. Harvey for 
their interest and encouragement during the course of this work. 



REFERENCES 

Benn|y, D. J. of Lin, C.C. 1960 On the secondary motion induced by oscil- 
lations in a shear flow, Phys. Fluids 4, 656. 

Hattia, F.R. & Nut ant, J. 1963 Detailed Flow field observations in the tran- 
sition process in^ a thick boundary layer. Proc. 1963 Heat Transfer and 
Fluid Mech . Institute , Stanford Univ. Press. — — 

Klebanoff , P.S. , Tidstrom, K.D. & Sargent, L.M. 1963 The three-dimensional 
nature of bound ary- layer instability. jJ. Fluid Mech . J2^ 1. 

Komoda, H. 1967 Nonlinear development of disturbance in a laninar boundary 
layer. Phys. Fluids . 10. SuppI.. 87. 

Mansour, N.N., Ferziger, J.H. & Reynolds, W.C. 1978 Large eddy simulation 
of a turbulent mixing layer. Dept. Mech. Engng . Stanford Uniy , Rep.TF-11. 

Milling, R.W. 1981 To 1 1 mi e n -Sc hi i cht i ng wave cancellation. Phys . Fluids . 

Fluid Mech J ^J 98 ^ 41 Nunerical investigation of turbulent channel flow, 

Moin, P., Reynolds, W.C. & Ferziger* J.H. 1978 Large eddy simulation of i n- 
compressible turbulent channel flow. De£t. Mech. Engng! Stanford Uniy 

Nishioka, M. , Asai, M. & Iida, S. 1980 An experimental investigation of the 
secondary instability. In Laminar-Turb ulent Transition (ed.R. EddI er of 
H. Fasel), Academic Press, T77 ~ 

Nishioka, M. , Asai , M. & Iida, S. 1981 Wall phenomena in the final stage of 
Acad^ni*c°Press TU 113 1 6nCe * In Transition and Turbulence (ed. R.E. Meyer), 

ln»l I »TC °L P n U n^ pecl, ' >) and ^ ct " al Station. 

Orszag, S.A. & Kells, L.C. 1980 Transition to turbulence in plane Poise- 
ui lie and plane Couette flow. J. Fluid Mech ., 96, 159. 

Orszag, S.A. & Patera, A.T. 1981 Subcritical transition to turbulence in 
planar shear flows. Phy . Rev , Let . 45 , 989. 

Reynolds, W.C. 1967 A Fortran IV program for solution of the Orr-^Sontmerfeld 
equation. Dept . Mech . Ehghg . Stanford Univ . Rep, fM4. 

Smith, C. R. & Metzleh, S. P. 1983 the characteristics of low- speed streaks 
129 th 27 near ~ Wa1 r69ioh of a turbulent boundary layer. J, Fluid Mech .. 


26 



Tani, I. 1969 Boundary layer transition. Ann. Rev. Fluid Mech . JU 169. 

Tennekes, H. & Lunley, J.L. 1972 A First Course in Turbulence . MIT Press. 

Theodorsen, T. 1955 The structure of turbulence in 50 Jahre Grenzschichtfo- 
rschung . Friedr. Vieweg of Sohn, 55. " 


27 








Amplitude 



Figure 2. Time-history of maximum disturbance amplitudes. 

(a) Two-dimensional primary; (b) three-dimensional 
primary; (c) two-dimensional harmonic. 
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Figure 9. Spanwiise variations Of u, . (a) T =* 0; (b) T = 17; 
(d) T *• 44; (e) T ** 79. (i) bottom, xp = -0.968; (ii) middle 

(iii) top, Xo ■« -0.764. * 
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Figure 12. Contour plots of a'uj./'Sx^. (a) Three-spike stage, figure 5 from 

Nishioka et al. (1981)'; (b)‘ computations at 1? * 27 in the x 1 -x 2 
plane, contours from 0.3 to 6.8. 





Figure 13. Contour plots of 3^/ax (a) five-spike stage, figut 

. Nishioka et al. f 1981) ; (b) computations at T =l 

in the x r x 2 plane; contours from 0.7 to 7.0. At thi< 
spanwise position, we define xorCxo) 
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